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Distinct population structure for co-occurring Anopheles goeldii 
and Anopheles triannulatus in Amazonian Brazil 
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To evaluate whether environmental heterogeneity contributes to the genetic heterogeneity in Anopheles triannu- 
latus, larval habitat characteristics across the Brazilian states of Roraima and Para and genetic sequences were ex- 
amined. A comparison with Anopheles goeldii was utilised to determine whether high genetic diversity was unique to 
An. triannulatus. Student t test and analysis of variance found no differences in habitat characteristics between the 
species. Analysis of population structure of An. triannulatus and An. goeldii revealed distinct demographic histories 
in a largely overlapping geographic range. Cytochrome oxidase I sequence parsimony networks found geographic 
clustering for both species; however nuclear marker networks depicted An. triannulatus with a more complex his- 
tory of fragmentation, secondary contact and recent divergence. Evidence of Pleistocene expansions suggests both 
species are more likely to be genetically structured by geographic and ecological barriers than demography. We 
hypothesise that niche partitioning is a driving force for diversity, particularly in An. triannulatus. 
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Malaria remains one of the most important tropical 
diseases in the world with an estimated 216 million cases 
in 2010 (WHO 2011). Brazil has the largest incidence 
of disease and malaria related deaths in South America, 
averaging 300 thousand cases annually (MS 2008, Bar- 
reto et al. 2009). Nearly all cases occur in the Amazon 
Basin (Akhavan et al. 1999, Moutinho et al. 2011), where 
the relatively low levels of transmission have been diffi- 
cult to control (da Silva-Nunes et al. 2011). The Amazon 
Basin is one of the world's most important bioregions, 
with a variety of habitats suitable fox Anopheles species, 
many of which are malaria vectors. 

Most primary vectors are widespread because of 
excellent colonising abilities and adaptation to variable 
environmental conditions (Costantini et al. 2009, Cohuet 
et al. 2010, Loaiza et al. 2012). Though not major ma- 
laria vectors in South America, Anopheles goeldii and 
Anopheles triannulatus are locally and regionally im- 
portant (Rubio-Palis 1994, Povoa et al. 2001, Galardo et 
al. 2007) with large geographic distributions and a mod- 
erate to wide habitat range (unpublished observations). 
Both species can colonise and increase abundance in 
altered or temporary environments (Tadei & Dutary- 
Thatcher 2000). Overlapping species distributions and 
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mostly generalist behaviour (McKeon et al. 2013), in ad- 
dition to the hypothesis of competition towards exclu- 
sion, noted by Tadei et al. (1998) make An. goeldii and 
An. triannulatus ideal comparisons. 

An. goeldii (Rozeboom & Gabaldon), often mistak- 
enly identified as Anopheles nuneztovari (Sallum et al. 
2008), has been resurrected from synonymy (Calado et 
al. 2008) and is part of the Nuneztovari complex (Mon- 
toya-Lerma et al. 2011). An. goeldii, thought to be re- 
stricted to the Amazon Basin (Calado et al. 2008), may 
extend further north in sympatry with An. nuneztovari 
(Scarpassa & Conn 2011). The close genetic association 
between An. nuneztovari, an important malaria vector in 
Colombia and Venezuela (Bourke et al. 2010), and An. 
goeldii, raises concerns over the putative vector status of 
the latter and its habitat specificity. Thus, studies involv- 
ing the population structure and ecological distribution 
of An. goeldii are useful in determining the extent of its 
distribution and divergence. 

An. triannulatus (Neiva & Pinto) has been reported 
from Nicaragua to Argentina (Faran 1980). Originally 
described as a polymorphic taxon (Faran & Linthicum 
1981), the intraspecific variation was attributed to ad- 
aptation to different habitats (Forattini 1962, Faran 
1980, Rosa-Freitas et al. 1998). More recently, it has 
been recognised as a complex of at least three species 
(Silva-do-Nascimento & Lourenco-de-Oliveira 2002, 
Silva-do-Nascimento et al. 2006, 2011). Recent stud- 
ies showing anthropophilic feeding and biting activity 
both indoors and out (Brochero et al. 2006, Gutierrez 
et al. 2009), support vector incrimination in some areas 
of Latin America (Silva-do-Nascimento et al. 2006, de 
Barros et al. 2007). 
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As factors affecting speciation and population dif- 
ferentiation evolve at different rates (de Queiroz 2007), 
local and regional analysis of ecology and population 
structure may identify behavioural or ecological dif- 
ferences that complement divergence studies. Use of 
nuclear markers in addition to mitochondrial markers 
reduces the process error associated with coalescence 
(Zink & Barrowclough 2008) and has led to an infor- 
mal set of genes sequenced frequently for phylogenetic 
studies including the mitochondrial cytochrome oxidase 
I (COI) (Avise 2000, Molina-Cruz et al. 2004, Sunghoon 
et al. 2011), nuclear white gene (Besansky & Fahey 1997, 
Reidenbach et al. 2009) and internal transcribed spacer 
(ITS)2 gene family (Zapata et al. 2007, Shultz & Wolf 
2009, Wiemers et al. 2009). Within the mitochondrial ge- 
nome the more conservative "folmer region" has become 
the standard barcode fragment with a 3% divergence 
threshold for speciation events (Hebert et al. 2003). 

The diversity of Neotropical anophelines has been 
attributed to the ability of dipterans to readily adapt to 
and utilise a broad variety of ecological niches (Grimaldi 
& Engel 2005) with more generalised populations being 
more ecologically heterogeneous (Bolnick et al. 2007). 
Thus, the genetic heterogeneity seen in the Triannulatus 
and Nuneztovari complexes may be explained, in part, 
by environmental heterogeneity (Gram & Sork 2001), 
including local or microgeographic adaptation, past and 
present ecological barriers and/or demographic events 
(Brouat et al. 2004). Additionally, information on the 
dynamics of larval anophelines is limited (Shililu et al. 
2003) and may provide insights into population differ- 
ences, thereby complementing adult studies. 

The purpose of this study was to examine the distri- 
bution, abundance and genetics of largely generalist and 
sympatric species from 76 larval sites to address the fol- 
lowing questions: (i) what is the level of habitat differen- 
tiation between species, (ii) do patterns of diversity in^4«. 
goeldii and An. triannulatus indicate a common cause 
and (iii) can genetic differentiation be explained by de- 
mographic phenomena and (or) geographic boundaries? 

MATERIALS AND METHODS 

Mosquito collection - Fourth instar anopheline lar- 
vae were collected from 13 localities, including eight in 
the state of Roraima (RR) in northern Amazonian Bra- 
zil and five in the state of Para (PA) in eastern Amazo- 
nian Brazil (Fig. 1, Supplementary data). Sampling was 
done once for each randomly selected larval site. If the 
site was positive for anophelines, collections continued 
for an hour to measure relative species abundance. The 
study was conducted over three years (2009-2011) dur- 
ing malaria transmission season in each state. The larval 
sites in RR were located between Amajari and Ecuador, 
roughly 500 km along the RR-174, in 2009 and 2011, and 
in PA between Moju and Maraba, approximately 560 km 
along PA-475/263/150 in 2010. The collection protocol 
was approved by the Evandro Chagas Institute, Para State 
Ethical Committee and the New York State Department 
of Health, Institutional Review Board. Physical charac- 
teristics of the larval habitats were recorded. These in- 
cluded biotic factors such as temperature, alkalinity, pH, 




Fig. 1: map of Brazil indicating the location of the 13 localities and 
species distributions. See Supplementary data for more information. 
PA: state of Para; RR: state of Roraima. 



conductivity, salinity and turbidity. Additionally, canopy 
coverage was estimated by spherical densitometer (Mi- 
nakawa et al. 2005). 

Larvae were identified morphologically as An. trian- 
nulatus s.l. and An. nuneztovari using the key of Deane et 
al. (1946) and confirmed by species specific ITS2 restric- 
tion fragment length polymorphism (Zapata et al. 2007). 
Morphological keys do not reflect the recent differences 
between the Colombian and Venezuelan An. nuneztovari 
and the Brazilian An. goeldii, samples, therefore samples 
were sequenced and compared to those from Calado et 
al. (2008) using a Bayesian inference (BI) phylogenetic 
approach (data not shown). Total genomic DNA was ex- 
tracted using the DNeasy Blue Tissue Kit (Qiagen, CA, 
USA) and maintained at Griffin Laboratory at -80°C. 

Statistical analysis - A paired sample t test was used 
to compare An. triannulatus and An. goeldii larval oc- 
currence across all sites. Subsequent analysis of simi- 
lar larval sites based upon the habitat type and type of 
water body sampled (pond, pool, swamp etc.), examined 
variation in larval densities using two-way analysis of 
variance (ANOVA). Additional analyses examining 
differences in habitat ecology, based on biotic factors, 
surrounding vegetation and canopy coverage were ex- 
amined to test the hypothesis that environmental hetero- 
geneity is positively correlated with genetic diversity. 

Amplification and sequencing - A representative 
sample of one -50 mosquitoes per locality was selected 
for DNA amplification and sequencing of the COI. Sub- 
sets of up to 10 samples per locality for the white gene 
and up to five for ITS2 were amplified and sequenced. 
The 65 8 -bp barcode fragment of the COI gene was am- 
plified using the forward primer LCO1490 and the re- 
verse primer HC02198 (Folmer et al. 1994). Individual 
polymerase chain reaction (PCR) reactions were per- 
formed using Ready-To-Go RT-PCR Beads (Amersham 
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Pharmacia/Biotech, NJ, USA) or PCR Supermix High 
Fidelity (GoTaq, Promega, CA) run on a Bio-Rad PTC- 
100 or 200 series thermal cycler (Bio-Rad Inc), using 
the conditions stipulated in Hebert et al. (2003). The Ap- 
plied Genomics Technology Core (Wadsworth Center) 
carried out the sequencing on an ABI PRISM 3700 au- 
tomated DNA sequencer. The forward and reverse COI 
sequences were aligned using Sequencher 3.0 (Gene 
Codes Corps, MI, USA), grouped together by sight and 
trimmed in MEGA v. 3.1 (Kumar et al. 2004) creating a 
600-bp fragment. In addition, a subset of three-five An. 
tricumulatus samples from each locality was sequenced 
following the protocol of Pedro and Sallum (2009), to 
yield a 449-bp fragment of the 3' COI. High levels of 
intraspecific variation, deep structure and complex phy- 
togenies based on the 3' end of the COI gene (Moreno et 
al. 2013) are not congruent with recently described spe- 
cies An. triannulatus, Anopheles haplophylus and An. 
triannulatus "C" (Silva-do-Nascimento & Lourenco- 
de-Oliveira 2002, Silva-do-Nascimento et al. 2006, 
2011). The barcode region is more variable than the 3' 
end of the COI (Ruiz-Lopez et al. 2012) and able to im- 
prove resolution at deeper nodes (Saunders & Le Gall 
2010). Therefore, analyses herein focused on the barcode 
fragment and included only a subset based on the 3' end 
to allow systematic comparisons of An. triannulatus 
populations across Brazil. 

Nuclear genes were included to investigate the his- 
tory of the nuclear genome and provide independent 
evidence of population divergence (Loiaza et al. 2012). 
Additionally, nuclear genes do not suffer the short com- 
ings of mitochondrial markers (haploid, small genome 
size, short coalescent time, maternal inheritance) (Bal- 
lard & Whitlock 2004) and are not as biased in base 
composition (Brower & De Salle 1998). An 800-bp frag- 
ment of the white gene was amplified using the W2R 
and WF primers, with the PCR conditions as reported 
in Mirabello and Conn (2008). A 500-bp fragment of 
the ribosomal ITS2 was amplified using the primers 18S 
and 28 S following the parameters in Li and Wilkerson 
(2005). The PCR products were cleaned, sequenced and 
aligned creating a 570 -bp fragment of the white gene in- 
cluding the intron. The final product size for ITS2 was 
433 bp and 371 bp for An. triannulatus and An. goel- 
dii, respectively. A single sequence per individual/gene 
was included in all subsequent analyses. Haplotypes of 
heterozygous individuals were inferred for each set of 
gene sequences using PHASE 2.1 software (Stephens et 
al. 2001, Stephens & Donnelly 2003, White et al. 2009). 
Unique sequences for all markers are available in Gen- 
Bank (accessions KC167670-KC167816). 

Genealogical relationship - Genealogical trees for 
each species were estimated for the concatenated (COI 
+ white) data set using BI analysis performed with Mr- 
Bayes v. 3.1 (Huelsenbeck & Ronquist 2001, Ronquist & 
Huelsenbeck 2003). Data were partitioned by gene using 
the model of nucleotide substitution (HKY+I and TrN+G 
for An. triannulatus and TrN+I+G and TPM2uf+G for 
An. goeldii) that best fit the white gene and COI, respec- 
tively, determined by jModelTest (Posada 2008). ITS2 
sequences were excluded from these analyses because of 



the limited sample size. The settings were two simultane- 
ous, independent runs of the Markov Chain Monte Carlo 
for two million generations, sampling every 1,000 gen- 
erations with a "burnin" of 25%. The outgroup Anophe- 
les marajoara (DQ076221 and AY956296 deposited in 
GenBank) was chosen based upon its phylogenetic posi- 
tion in Sallum et al. (2000). 

Genetic variation - Genetic structure of multiple 
populations based on statistical parsimony groupings 
were examined by analysis of molecular variance (AM- 
OVA) using Arlequin v. 3.1.1 (Excoffier et al. 2005) and 
spatial analysis of molecular variance (SAMOVA), v. 1.0 
(Dupanloup et al. 2002). Both algorithms were used to 
cluster the 600-bp barcode fragment data into geneti- 
cally and in the case of SAMOVA, geographically, ho- 
mogeneous populations by generating F statistics (F sc , 
F sp F CJ ,). Each species was analysed separately, An. 
triannulatus into groups of K = 2-16 and An. goeldii 
into groups of K — 2-20, with 1,000 simulated anneal- 
ing steps from each of 100 sets of initial starting condi- 
tions. Groups were based on individual larval sites. A 
second analysis was undertaken for An. triannulatus us- 
ing the 449-bp fragment haplotypes of Pedro and Sallum 
(2009) (GU445849-899 and AF417702) combined with 
the samples in the present study, to test their hypothesis 
of geographical sub-structuring within the species based 
on five populations across Brazil. 

Population structure and demographic history - A 
statistical parsimony network (SPN) estimated genealog- 
ical relationships among haplotypes with a 95% identity 
for each marker using TCS 1.13 (Clement et al. 2000), 
except An. goeldii ITS2 sequences, which were exam- 
ined at 94% sequence similarity. Homoplasy in all net- 
works was resolved using the algorithm estimation rules 
in Crandall and Templeton (1993). Estimates of time to 
coalescence were calculated for the COI fragment only 
and compared using 6 S values (Watterson 1975). 

The differentiation and polymorphism statistics for 
COI and white gene sequences by species and local- 
ity were computed in DnaSP, v. 5.0 (Librado & Rozas 
2009). Because of the high proportion of shared haplo- 
types among habitats within a single locality, samples 
from an individual town/city were pooled to more accu- 
rately represent diversity statistics. Furthermore, haplo- 
type distributions across four ecoregions were examined 
for each of the seven environmental variables to test the 
hypothesis that climatic differences sustain (and drive) 
diversity (Hoorn et al. 2010). To test this hypothesis, we 
considered both the availability of breeding sites and 
haplotype diversity as measures of gene flow. Individual 
sites were grouped within ecoregions, based on descrip- 
tions from Rosa-Freitas et al. (2007). 

The hypothesis of strict neutrality for each species 
was examined using DnaSP, with the statistics D T (Ta- 
jima 1989) and D and F (Fu 1993) evaluating the likeli- 
hood of background selection and Fu's F (Fu 1997) and 
R 2 (Ramos-Onsins & Rozas 2002) examining evidence 
of population expansion. The mismatch distributions 
(simulated in Arlequin) were conducted to confirm de- 
mographic expansion using sudden and spatial models of 
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expansion. Statistically significant differences between 
observed and simulated distributions were evaluated 
with the sum of square deviations to reject the hypoth- 
esis of demographic expansion (Zarza et al. 2008). Dates 
of population expansion were estimated for the COIgem 
with the formula T= x/2u (Rogers & Harpending 1992) 
using 10 generations per year and the mutation rate of 
Drosophila estimated at 1 x 10 8 substitutions per site 
per year (Walton et al. 2000, Yang et al. 2011). 

RESULTS 

Site diversity and larval abundance - Larvae were 
predominantly collected from the edges of naturally oc- 
curring ponds, comprising 49.19% (n = 200) of the total 
Anopheles identified. An. goeldii and An. triannulatus 
larvae were relatively uncommon in collections from 
seepages, stream margins and swamps, where they com- 
prised only 10.9% of the total. Both species were col- 
lected from the same site 17 times (Supplementary data). 
Twelve sites contained An. triannulatus and not An. 
goeldii, whereas 13 sites contained An. goeldii, but not 
An. triannulatus. A large proportion, 68.97% (n = 140), 
of An. triannulatus larvae were sampled from ponds. 
An. goeldii larvae appeared equally likely to inhabit ar- 
tificial pools (37.93%, n = 77), ponds (29.55%, n = 60) 
and temporary ditches (25.12%, n = 51). No correlations 
between abundance and habitat type for either species 
were detected by ANOVA (data not shown). 

Habitat characteristics of An. triannulatus and An. 
goeldii were similar: the 76 sites were not differentiated 
by larval abundance (t = -0.2980, df= 42, p = 0.3836). 
ANOVA detected a significant difference in turbid- 
ity between sites containing only An. goeldii and those 
negative for either species (F = 4.89, p = 0.032) that was 
not supported by a t test (t = 1.179, df = 16, p = 0.055). 
There are no other significant differences between en- 
vironmental variables for sites unique for either species 
and sites either positive or negative for both (Table I). 

Genealogical relationship - The concatenated (COI 
+ white gene with intron) BI tree for An. triannulatus, 
with 50 parsimoniously informative sites, indicated two 
strongly supported (97 and 98 posterior probabilities) 



distinct populations of An. triannulatus (Fig. 2). The 
two populations correspond to the Brazilian states, with 
the exception of samples from Moju and Goianesia (9 
and 10, respectively), that were included in the northern 
group (RR localities). The close relationship between 
the RR samples and localities Moju and Goianesia was 
also evident in the SPN. This may indicate a historical 
boundary separating north and south populations, with 
the northern limit occurring just south of the Amazon 
River and including parts of Goianesia and Moju. The 
more derived population 1 corresponds to localities in 
the Amazon Basin, near the hypothesised population 
origin (Pedro & Sallum 2009), with the majority of pop- 




Fig. 2: concatenated (cytochrome oxidase I + white gene) Bayesian 
inference trees, partitioned by gene using the model of nucleotide 
substitution that best fit the data as determined by jModelTest. blue: 
samples from Goianesia, Tucurui, Jacunda and Maraba [state of Para 
(PA)] corresponding to An. triannulatus population 3; green: Moju 
and Goianesia samples (PA) in addition to An. triannulatus popula- 
tion 2; orange: localities in the state of Roraima and An. triannulatus 
population 1. Asterisk means Anopheles marajoara outgroups. 



TABLE I 

Characteristics (± 95% confidence intervals) of habitats with Anopheles triannulatus only, Anopheles goeldii only, 

both species or neither species present 





An. triannulatus 


An. goeldii 


Both 


Neither 


Collections (n) 


12 


13 


17 


34 


pH 


5.5 ±0.00 


5.58 ±0.15 


5.68 ±0.24 


5.62 ±0.16 


Alkalinity (ppm) 


35 ± 29.40 


33.85 ± 27.50 


47.06 ± 27.91 


52.35 ±42.58 


Temperature (°C) 


26.29 ± 1.37 


27.44 ± 1.88 


26.79 ± 1.10 


26.48 ± 0.95 


Conductivity (|j,S/cm) 


35.68 ± 16.33 


39.49 ± 19.68 


34.01 ± 12.93 


26.04 ± 9.76 


Salinity (ppm) 


18.09 ±8.13 


19.92 ± 10.31 


25.11 ± 17.34 


13.24 ±4.90 


Turbidity (JTUs) 


21.67 ± 15.93 


28.46 ± 16.88 


17.50 ±6.88 


13.09 ±5.36 


Canopy coverage (%) 


11.28 ± 14.20 


13.74 ± 13.43 


13.77 ± 11.38 


10.55 ±6.33 
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ulation 2 occurring along one of the proposed routes of 
lineage expansion. Within grouping 2, smaller subdivi- 
sions, related to geography and corresponding to popula- 
tions 2 and 3, were apparent between haplotypes 9 and 3 
that are separated by six mutations and two haplotypes 4 
and 1, separated by two mutations (Fig. 3). 

The concatenated (COI+ white gene with intron) BI 
tree for An. goeldii, with 32 parsimoniously informa- 
tive sites, indicated two moderately supported (72 and 
68 posterior probabilities) populations which correspond 
to the Brazilian states (Fig. 2). Unlike An. triannulatus, 
there does not appear to be any sub-structuring in the 
samples from PA or support for a separate Goianesia/ 
Mojii population. 

Genetic variation - AMOVA indicated consistently 
higher F ST and F CT fixation indices across genes and 
various COI groupings for An. triannulatus, i.e., ge- 
netic structuring among and within localities rather than 
between groups (Table II). However, there was strong 
genetic differentiation (41.01%) also seen in the COI 
grouping between states, which correlates well with the 
groups defined by the Bayesian tree. 

The small genetic differentiation of An. goeldii be- 
tween localities suggests panmixia with no significant 
barriers to gene flow within individual states. AMOVA 
indicated a large proportion of variation between the two 
states and among samples from different breeding site 
types in a given locality (Table II). Genetic structuring 
in the latter may indicate environmental differences and 
niche partitioning as drivers of diversity. 



SAMOVA analysis for both species found no evidence 
of a geographic barrier. The addition of the An. triannula- 
tus 449-bp fragment of the COI gene (n = 35) to all sam- 
ples from Pedro and Sallum (2009) was unable to replicate 
the five groupings originally described across Brazil. 

Population structure and demographic history - 
Analysis of mitochondrial data from 129 An. triannu- 
latus and 133 An. goeldii specimens (Fig. 3) led to the 
discovery of two very different demographic histories 
between these sympatric species, which often co-occur 
within the same larval habitats. The time of the begin- 
ning of lineage divergence was calculated following 
T>J2k, where D a is the average nucleotide divergence 
between lineages and 2k is the divergence rate (Heckel 
et al. 2005). We utilised the mutation rate of Drosophila 
(Walton et al. 2000, Yang et al. 2011) to estimate these 
parameters and their significance based on the coa- 
lescence (DnaSP v. 5.10) (Librado & Rozas 2009) and 
MEGA v. 5.1 (Tamura et al. 2011). 

A parsimony network for An. triannulatus indicated 
that a minimum of six mutational steps separate three 
potential populations based on geography (Fig. 3). The 
most common haplotypes were 1 (n = 18), 2 (n = 15), 
3 (n = 11) and 6 (n = 9) (Fig. 3), with seven haplotypes 
(15.6%) shared among localities and 38 (84.4%) unique 
(Supplementary data). Overall, An. triannulatus was 
consistent with a category I phylogeographic pattern 
(Avise 2000) characterised by pronounced genetic gaps 
between spatially circumscribed haplogroups. A high 
proportion of singletons (23, 51.1%), in addition to in- 




Fig. 3: statistical parsimony network of cytochrome oxidase I sequences at 95% support, black: localities in the state of Roraima; grey: Moju and 
Goianesia samples from the state of Para (PA); n: number of sequences; white: samples from Goianesia, Tucurui, Jacunda and Maraba (PA). For 
An. triannulatus, black refers to population 1, grey (with the addition of haplotype 14) to population 2 and white to population 3. 
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TABLE II 

Analysis of molecular variance calculations for inter and intravariation among species and populations for various markers 



Anopheles triannulatus Anopheles goeldii 







Variance 


Fixation 


Variance 


Fixation 


Gene 


Source of variation 


(%) 


index 


(%) 


index 


COI 


Between RR and PA 


41.01 


F CT = 0.41" 


45.12 


F CT = 0.45" 




Among localities within states 


29.23 


F sc = 0.50" 


5.29 


F sc = 0.10* 




Within a given locality 


29.76 


F ST = 0.70* 


49.59 


F ST = 0.50* 


white 


Between RR and PA 


-3.01 


F CT = -0.03 1 - 


89.46 


F CT = 0.89" 




Among localities within states 


58.62 


F sc = 0.57* 


0.42 


F sc = 0.04' 




Within a given locality 


44.38 


F ST = 0.56" 


10.12 


F ST = 0.90* 


COI 


Between SPN populations 


38.85 


F CT = 0.39' 








Among localities 


33.92 


F sc = 0.55* 








Within a given locality 


27.23 


F ST = 0.73* 







a: p < 0.05; b: p < 0.001; c: not significant; COI: cytochrome oxidase I; PA: state of Para; RR: state of Roraima; SPN: statistical 
parsimony network derived groupings. 



dividual populations within An. triannulatus containing 
star-shaped nodes (surrounding haplotypes 2, 3 and 5), 
supports demographic expansions, background selec- 
tion or selective sweeps (Slatkin & Hudson 1991, Fu 
1997). Based upon the 600-bp fragment, D between the 
northern and southern states is 0.0118 [standard devia- 
tion (SD) = 0.003]. Therefore, the estimated divergence 
is 437,500-737,500 years ago, in the Pleistocene. The two 
PA populations had a D a = 0.014 (SD = 0.004), indicat- 
ing more ancestral divergence (congruent with Bayesian 
tree), approximately 500,000-900,000 years ago. 

An. goeldii COI is consistent with a category III phy- 
logeographic pattern (Avise 2000) where most haplotypes 
are closely related, yet some geographic localisation ex- 
ists. The most common haplotypes were 46 (n = 21), 47 (n 
= 33) and 48 (n = 13) (Fig. 3), with 14 haplotypes (41.2%) 
shared among localities, and 20 (58.8%) unique (Supple- 
mentary data). Similar to An. triannulatus, there was a 
relatively high proportion of singletons (14, 41.2%), in- 
dicating a demographic expansion, background selection 
or selective sweep (Slatkin & Hudson 1991, Fu 1997), 
particularly in PA (singletons arising from haplotype 56). 
In contrast, haplotypes from RR exhibited longer branch- 
es, missing haplotypes and a near equal distribution of 
shared haplotypes and single mutations, consistent with 
a signal of an older lineage with balancing selection (Fig. 
3). D between northern and southern states is 0.008 (SD 
= 0.002). Therefore, the estimated divergence is 316,000- 
516,500 years ago (Pleistocene). 

Networks of An. triannulatus nuclear data (white 
gene and ITS2) indicate a history of contractions and 
expansions. Similar to the COI gene, the ITS2 network 
depicts a category I pattern whereas the white gene re- 
veals a much more homogeneous haplotype distribution. 
Overall, the white gene follows a category II intraspecific 
pattern, with pronounced gaps between some branches 
and some clustering of anciently separated geographic 
populations (Fig. 4). 




An. triannulatus 

n = 25 

Fig. 4: statistical parsimony network of white gene (A) and internal 
transcribed spacer 2 sequences (B) at 95%, except ITS2 An. goeldii 
which is at 94% sequence support, black: localities in the state of Ro- 
raima; grey: Moju and Goianesia samples from the state of Para (PA); 
n: number of sequences; white: samples from Goianesia, Tucurui, Ja- 
cunda and Maraba (PA). 



Networks of An. goeldii nuclear data (white gene and 
ITS2) reveal distinct geographical populations. The white 
gene has a distinct category I pattern separated by five 
mutational steps, while the ITS2 haplotype clusters (Fig. 
4) are less pronounced, but still correspond primarily to 
RR and PA with a few outliers. Polymorphism analyses 
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of the CO/sequences found haplotype and nucleotide di- 
versity was greater among An. triannulatus populations 
(Supplementary data) in contrast to An. goeldii. 

Genetic polymorphism analyses of the COI sequenc- 
es indicated that haplotype and nucleotide diversity was 
greater among An. triannulatus populations (Supple- 
mentary data) in contrast to An. goeldii. The white gene 
exhibited lower nucleotide diversities compared to COI, 
but surprisingly An. triannulatus continued to demon- 
strate high haplotype diversity (Supplementary data), 
which combined with relatively low nucleotide diversity 
may signal a rapid expansion (Avise 2000, de Jong et al. 
2011). High COI variation between states and geograph- 
ic populations is strongly supported by population dif- 
ferentiation statistics (Table III), although the G ST values 
were low, suggesting a history of gene flow (Table III). 

In general, neutrality tests were negative for both 
species (Table IV), suggesting an excess of rare poly- 
morphisms consistent with either negative/positive se- 
lection or an increase in population size, although only 
one of these tests (F) was significant for An. triannu- 
latus, i.e., populations 2 and 3. Mismatch distribution 
for sudden expansions did not support this finding and, 
instead, revealed marginally significant expansion pat- 
terns for only population 1, corresponding to RR. Sud- 
den and spatial populations for An. goeldii could not be 
rejected (Fig. 5). Estimated expansions for both species 
are during the Pleistocene, An. triannulatus 129,693 
years ago (59,857-194,877) and An. goeldii 112,387 years 
ago (17,111-187,449). 

DISCUSSION 

Distribution patterns of mosquito species are depen- 
dent on the availability of suitable aquatic habitats (Grillet 
2000). Overlapping geographic distributions, in addition 
to the absence of significant difference in species abun- 
dance across sampled habitats, suggest that An. goeldii is 
a suitable comparison species for An. triannulatus. The 
lack of any significant differentiation between sympatric 
habitats, to those unique to a species or absent of either, 
suggests these species may be more broadly adapted 
(Eisenberg et al. 2000, Chase & Knight 2003). However, 
co-occurrence of species does not imply the same ability 
to exploit shared habitats (Diabate et al. 2005). 

The combined sequence results support distinct de- 
mographic histories for the two species despite similar 
distribution and larval habitat co-occurrence. Across 
molecular markers An. goeldii consistently depicts two 
moderately supported populations corresponding to the 
two Brazilian states, however, An. triannulatus appears 
to have undergone historical fragmentation, subsequent 
secondary contact and homogenisation followed by more 
recent divergence due to allopatry. Either environmental 
changes (e.g., habitat availability and climatic oscilla- 
tions) alone cannot account for population divergence or 
differences in biology allow these two species to respond 
differently to climatic oscillations. 

The clear separation of populations based on Brazil- 
ian states in An. goeldii strongly supports the alternative 
hypothesis of geographic isolation driving divergence 
rather than environmental changes. While often contest- 



ed in the neotropics, the Pleistocene refugia hypothesis 
(Haffer 1969) postulates that climatic fluctuations result- 
ed in the fragmentation of continuous forest into refuges 
separated by expanses of grass-dominated savannah or 
desert. Alternatively, the Miocene marine incursion hy- 
pothesis suggests sea level changes during the late Ter- 
tiary (23.7-1.8 million years ago) created three isolated 
regions: the Brazilian Shield, the Guiana Shield and the 
eastern Andean slopes (Webb 1995, Bates 2001, Hoorn 
et al. 2010). The last commonly explored mechanism of 
divergence is the riverine barrier hypothesis, where ma- 
jor rivers act as barriers to gene flow, thus promoting 
speciation between populations on either side (Wallace 
1852, Aleixo 2004, Mirabello & Conn 2008). 

Although no boundary was indicated by SAMOVA, 
this does not discount the possibility that an historical 
boundary existed, separating north and south popula- 
tions. Although one study found evidence for waterways 
over 4 km wide being complete barriers to gene flow 
(Pedro & Sallum 2009), others have determined that the 
Amazon is not an effective barrier (Fairley et al. 2002, 
Mirabello & Conn 2008) and that geographic distance 
and/or differing demographic histories may be the main 
forces responsible for partitioning genetic variation 
(Mirabello & Conn 2008). 

The discrepancy regarding the inclusion of Mojii and 
some of Goianesia among An. triannulatus RR popula- 
tion, but their inclusion with geographically similar PA 
samples in An. goeldii may be attributed to individual 
species' responses to dispersal opportunities (Brouat et 
al. 2004). The two An. triannulatus populations identi- 
fied in Goianesia could also be the result of a suture zone. 



TABLE III 

Inter and intrapopulation differentiation of 
Anopheles triannulatus and Anopheles goeldii populations 
based on statistical parsimony network 



Comparison 





An. triannulatus 


An. goeldii 






RR vs. PA 


Three populations 


RR vs. PA 


H s 


0.89" 


0.86" 


0.86" 


K* 


1.77" 


1.41" 


1.60" 


Z* 


7.43" 


7.00" 


7.83" 


s 

nn 


0.99" 


1.00" 


0.97" 


f 


129" 


258" 


125.84" 


G ST 


0.06 


0.09 


0.045 


K, 


10.09 


10.09 


7.43 



a: p < 0.0001; G ST : genetic differentiation; H s : genetic differen- 
tiation based on haplotype data; K s *: differentiation based on 
sequence data; K: average number of nucleotide differences; 
n: number of total individuals compared; PA: state of Para; 
RR: state of Roraima; S : measures how often the nearest 
neighbours of sequences are found in the same population; / : : 
genetic differentiation based on allele frequencies; Z*\ rank 
statistic to analyse sequence similarity. 
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TABLE IV 

Multiple statistical analyses evaluating the probability for a recent population expansion 



Species 


n 


D T 




pa 


F s 


R 2 


Anopheles triannulatus 


129 


-0.22187 


-1.28038 


-0.99732 


-10.329° 


0.08650 


Roraima 


57 


-0.37846 


-1.32424 


-1.17022 


-1.488 


0.09200 


Para 


72 


-0.47543 


-0.97094 


-0.93438 


-5.842 


0.08900 


Population 1 


56 


-0.08718 


-1.11433 


-0.88956 


-1.209 


0.10200 


Population 2 


36 


-1.06921 


-0.55581 


-0.85324 


-3.582° 


0.07770 


Population 3 


37 


-1.02686 


-0.40080 


-0.71980 


-4.809" 


0.07800 


Anopheles goeldii 


133 


-0.30013 


-1.24177 


-1.01960 


-5.867 


0.08260 


Roraima 


103 


0.06971 


-1.6295 


-1.1537 


-1.677 


0.0989 


Para 


30 


-0.0415 


0.15645 


0.11023 


-3.321 


0.1185 



a: significance, p < 0.50; D: Fu and Li's D (Fu 1993); D.-. Tajima's D (Tajima 1989); F: Fu and Li's F (Fu 1993); F s : Fu's Fs (Fu 1997); 
n: sample size; R 2 : Ramos-Onsin's and Roza's i?, (Ramos-Onsins & Rozas 2002). 



a c 




Differences (n) ■ Observed ■ Model 



Fig. 5: representative mismatch distributions with calculated % and 
estimated time of sudden expansion. Bars indicate observed values. 
Line represents the model for all Anopheles triannulatus (A), all 
Anopheles goeldii (B), An. triannulatus population 2 (C) (Mojii and 
Goianesia samples) and An. triannulatus population 3 (D) (Goianesia, 
Tucurui and Jacunda samples), ya: years ago. 



Species dispersal activity and differential selection can 
cause hybrid zones to appear mobile and settle as bound- 
aries between environments or in regions of low popula- 
tion density (Barton & Hewitt 1985, Mallet 1993, 2010). 
Suture zones that span several species suggest divergence 
was driven by a common influence of environmental 
change (Morgan et al. 2010). Further evaluation of spe- 
cies in the Goianesia area could resolve the contribution 
of environmental change to population diversity. 

Overlapping geographic distributions, in addition to 
the absence of significant difference in species abun- 
dance across sampled habitats, implies little habitat 
differentiation between species and a broad ecological 
range. Although biological differences may drive di- 
vergence uniquely in each of these sympatric species, a 
common cause for diversity patterns cannot be ruled out. 



Geographic boundaries/distance, variation in dispersal 
activity and ecological range may each contribute to the 
distinct patterns of divergence seen in An. triannulatus 
and An. goeldii. As species complexes are common in 
anophelines, additional studies are needed to resolve the 
extent of divergence. Continued studies using multiple 
markers at a local and continental scale will resolve the 
status of speciation and may clarify what role, if any, 
a given population, lineage or species has on malaria 
transmission, because population differences influence 
peak biting times, feeding preferences and vectorial ca- 
pacity (Lounibos & Conn 2000). 
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TABLE 

Species distribution of anopheline larvae in different larval habitats by locality over three transmission seasons (2009-2011) 

Anopheles triannulatus Anopheles goeldii 



Locality 


Habitat type 


(n) 


(n) 


Amajari (1) 


Marsh 


0 


0 




Pool 


1 


1 


Alto Alegre \l ) 


uitcn 


U 


U 




Marsh 


4 


U 




rono 


n 
U 


7 




Seepage 


0 


u 




Stream margin 


I 


1 


Boa Vista (3) 


Ditch 


0 


0 




Marsh 


u 


u 




Pond 


u 


1 

4 




Pool 


U 


1 A 

1 4 




Seepage 


u 


0 




Stream margin 


u 


A 

u 


Mucajai (4) 


Ditch 


u 


A 

u 




r 001 






Iracema (5) 


Ditch 


1U 


43 




r ona 




1 1 


Petrohna (6) 


Ditcn 


U 


A 

U 




Pond 


U 


A 
U 




Stream margin 


U 


1 

J 


Martins Pereira (7) 


Stream margin 


i 
1 


I 




Swamp 


14 


A 

4 


Ecuador (8) 


Pond 


U 


I 




Seepage 


U 






Stream margin 


n 
u 


A 

u 




Swamp 


U 


A 

u 


Moju (9) 


Pond 


19 


1 A 

1 9 




Marsh 


DO 


5 1 




Stream margin 


U 


A 
U 


Goianesia (10) 


Marsh 


0 


A 

0 




Pond 


14 


A 
U 




Pool 


0 


A 

u 




Stream margin 


U 


A 

u 




Swamp 


0 


0 


Tucurui (11) 


Pond 


7 


2 




Stream margin 


6 


0 


Jacunda (12) 


Marsh 


5 


0 




Pond 


11 


0 




Stream margin 


0 


2 


Maraba (13) 


Ditch 


0 


6 




Marsh 


0 


0 




Pond 


2 


3 




Stream margin 


1 


0 


Total 




203 


230 



1-13: see Fig. 1 for more information. 
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TABLE 



Description of shared cytochrome oxidase I haplotypes and genetic polymorphism statistics 
for Anopheles triannulatus and Anopheles goeldii 





Site 


n 


Haplotypes 


;i (SD) 


1-Yf 2 (SD) 


P 


K 


An. triannulatus 


1 


1 


1 


- 


- 


- 


- 




2 


4 


1 (2), 40, 41 


0.01 (0.00) 


0.83 (0.22) 


12 


7.67 




5 


37 


1 (15), 2 (6), 7, 8, 9 (2), 35, 36 (2), 37, 39 (4), 43, 45 


0.01 (0.00) 


0.80 (0.05) 


20 


5.61 




7 


15 


2 (9), 14, 38, 42 (3), 44 


0.00 (0.00) 


0.63 (0.13) 


14 


2.27 




9 


30 


3 (11), 4 (7), 10, 11, 30, 31 (2), 32, 33, 34 (4) 


0.00 (0.00) 


0.81 (0.05) 


12 


1.95 




10 


14 


6 (3), 12, 13 (3), 15, 16, 17 (2), 24, 25, 28 


0.01 (0.00) 


0.92 (0.05) 


20 


6.99 




11 


13 


5 (2), 6, 16, 18, 19, 20, 21, 22, 26, 27 


0.01 (0.00) 


0.96 (0.04) 


16 


4 




12 


15 


5 (4), 6 (4), 23 (2), 24, 25 (2), 29 


0.01 (0.00) 


0.83 (0.06) 


9 


2.99 






129 




0.02 (0.00) 


0.95 (0.01) 


57 


10.1 


An. goeldii 


2 


4 


46, 47 (3) 


0.00 (0.00) 


0.50 (0.27) 


2 


1 




3 


15 


46 (4), 47 (7), 48, 49, 77, 78 


0.01 (0.00) 


0.74 (0.09) 


15 


3.26 




4 


18 


46 (3), 47 (6), 48 (2), 49 (2), 50, 63, 67, 72 


0.01 (0.00) 


0.86 (0.06) 


18 


6.78 




5 


50 


46 (7), 47 (14), 48 (10), 49 (3), 50 (2), 51, 52, 63 (3), 
70, 71,73(2), 74(2), 75 


0.01 (0.01) 


0.86 (0.03) 


24 


6.29 




6 


3 


47 (2), 72 


0.01 (0.01) 


0.67 (0.31) 


11 


7.33 




7 


6 


46, 50, 68, 69, 73, 76 


0.01 (0.00) 


1.00 (0.10) 


18 


8.47 




8 


7 


46 (5), 47, 53 


0.01 (0.00) 


0.52 (0.21) 


12 


3.81 




9 


17 


53, 56, 57, 59, 60 (2), 61 (3), 62 (2), 64, 65, 66 


0.00 (0.00) 


0.93 (0.04) 


10 


2.18 




11 


2 


54, 63 












13 


11 


54 (3), 55 (2), 58, 60, 63 (2), 64 


0.01 (0.00) 


0.91 (0.07) 


13 


5.6 






133 




0.01 (0.01) 


0.90 (0.002) 


44 


7.43 



K: average number of differences; n: number of sequences; P: polymorphic sites; SD: standard deviation; 1-Zf, 2 : haplotype di- 
versity; IT: nucleotide diversity; 1-13: see Fig. 1 for more information. Bold indicates shared haplotypes and indicates diversity 
among lineages or species. 



TABLE 



Description of shared white haplotypes and genetic polymorphism statistics for Anopheles triannulatus and Anopheles goeldii 



Species 


Site 


n 


Haplotypes 


7t(SD) 


;-I/f(SD) 


P 


K 


An. triannulatus 


1 


1 


B 












2 


3 


B, E (2) 












5 


10 


B (3), C (4), D (2), E 


0.00 (0.00) 


0.51 (0.16) 


1 


2.47 




7 


4 


E (2), P, T 












9 


9 


A (8), U 


0.00 (0.00) 


0.22 (0.17) 


1 


0.22 




10 


9 


E, N (2), O, Q (3), R, S 


0.01 (0.00) 


0.81(0.12) 


9 


3.94 




11 


10 


E, F, I, J (2), K, L (2), N (2) 


0.01 (0.00) 


0.93 (0.06) 


11 


4.02 




12 


8 


G (2), H, I, J, L (2), M 


0.00 (0.00) 


0.93(0.08) 


4 


1.54 






54 




0.01 (0.00) 


0.93 (0.02) 


24 


5.43 


An. goeldii 


2 


4 


V, W (2), CC 












3 


10 


V (3), W (6), X 


0.00 (0.00) 


0.600(0.13) 


4 


1.07 




4 


10 


V, W (9) 


0.00 (0.00) 


0.20 (0.15) 


1 


0.20 




5 


9 


W(9) 


0.00 (0.00) 


0.000(0.00) 


0 


0 




8 


1 


W 












9 


10 


Y (9), Z 


0.00 (0.00) 


0.200(0.15) 


2 


0.4 




13 


8 


Y (5), AA (2), BB 


0.00 (0.00) 


0.61 (0.16) 


3 


1.04 






52 




0.01 (0.00) 


0.66 (0.05) 


12 


2.85 



K: average number of differences; n: number of sequences; P: polymorphic sites; SD: standard deviation; 1-Sf 2 ; haplotype diver- 
sity; n: nucleotide diversity; 1-13: see Fig. 1 for more information. Bold indicates diversity among lineages or species. 



